﻿#pragma once
#include <pcl/filters/filter.h>

#include <memory>
#include <string>
#include <vector>
#include <utility>

using std::string;
using std::vector;
using std::pair;

namespace oil
{
class RasterSurface;

/**\brief ASF 滤波算法阈值
*
* \code 
* AdaptiveSurfaceFilter<PointXYZRGBL> surface_filter;
* ASFParameter param;
* param.initUrban();

* surface_filter.setInputCloud(cloud);
* surface_filter.setParameters(param);
* surface_filter.filter(*ground);
* \endcode
*/
struct ASFParameter
{
	///\brief 每级金字塔窗口大小系数，乘以step_size
	double step_size_;

	///\brief 最后阈值t = bend_gain * bend_mask +scale_gain + threshold
	double threshold_;

	///\brief 最粗的金字塔的窗口大小
	double max_window_;

	///\brief 最精细的金字塔的窗口大小
	double min_window_;

	///\brief 最大的scale gain，从粗到精会依次减小
	double max_scale_gain_;

	///\brief 最大的bending gain，根据bending energy变化
	double max_bending_gain_;

	///\brief 最大的正则化系数，从粗到精会增加
	double max_lambda_;

	///\brief 分段内插的系数
	int bending_level_;

	///\brief radius outlier filter 需要设置窗口大小，为min window * m
	int m_;

	///\brief radius outlier filter 需要设置窗口中的点数，少于k的认为是outlier
	int k_;

	///\brief 是否检测low point
	bool detect_low_outlier_;

	///\brief low outlier detection的参数
	///\sa LowOutlierFilter
	double low_outlier_mophological_window_;
	double low_outlier_threshold_;
	double low_outlier_mophological_size_;
	int max_low_grids_;

	///\brief 在建筑物不是特别高的城区，这个地方low outlier 相对少，且不集中，并且能适当处理一些地形变化
	void initUrban()
	{
		m_ = 5;
		k_ = 3;
		max_window_ = 60;
		threshold_ = 0.3;
		max_bending_gain_ = 0.0;
		min_window_ = 2.0;
		max_scale_gain_ = 0.3;
		step_size_ = 1.4;
		max_lambda_ = 0.5;
		detect_low_outlier_ = true;

		//相当于是25m的范围内的中值
		low_outlier_mophological_size_ = 2;
		low_outlier_mophological_window_ = 5.0;
		low_outlier_threshold_ = 0.5;
		max_low_grids_ = 6; //这个现在被中值替代掉
	}

	///\brief 在山区（植被），山区通常很少有low outlier
	void initMountain()
	{
		m_ = 5;
		k_ = 3;
		max_window_ = 30;
		threshold_ = 0.3;
		max_bending_gain_ = 0.3;
		min_window_ = 1.0;
		max_scale_gain_ = 0.3;
		step_size_ = 1.4;
		max_lambda_ = 0.5;
		detect_low_outlier_ = false;
	}

	///\brief 在建筑物较大，较高的区域，通常low outlier 很多，而且密集，但是这些区域通常不太复杂，地形平坦
	void initMetropolitan()
	{
		m_ = 5;
		k_ = 3;
		max_window_ = 80;
		threshold_ = 0.3;
		max_bending_gain_ = 0.0;
		min_window_ = 2.0;
		max_scale_gain_ = 0.3;
		step_size_ = 1.4;
		max_lambda_ = 0.5;
		detect_low_outlier_ = true;

		//即四周13*5的范围内的中值
		low_outlier_mophological_window_ = 5.0;
		low_outlier_threshold_ = 0.5;
		low_outlier_mophological_size_ = 6;
		max_low_grids_ = 6;
	}
	ASFParameter(){
		step_size_ = 1.2;
		threshold_ = 0.3;
		max_window_ = 30.0;
		min_window_ = 1.0;

		max_scale_gain_ = 0.3;
		max_bending_gain_ = 0.5;
		max_lambda_ = 0.5;

		bending_level_ = 5;

		m_ = 5;
		k_ = 3;

		detect_low_outlier_ = false;

		low_outlier_mophological_window_ = 2.0;
		low_outlier_threshold_ = 0.5;
		low_outlier_mophological_size_ = 2;
		max_low_grids_ = 6;
	}
};

/**\brief ASF进行滤波
*
*采用金字塔影像策略，进行滤波，详见：
Hu, H., Ding, Y., Zhu, Q., Wu, B., Lin, H., Du, Z., Zhang, Y. and Zhang, Y., 2014. An adaptive surface filter for airborne laser scanning point clouds by means of regularization and bending energy. ISPRS Journal of Photogrammetry and Remote Sensing, 92: 98-111.
*/
template <typename PointT>
class PCL_EXPORTS AdaptiveSurfaceFilter : public pcl::Filter < PointT >
{
protected:
	using Filter<PointT>::filter_name_;
	using Filter<PointT>::getClassName;
	using Filter<PointT>::input_;
	using Filter<PointT>::indices_;

	typedef typename Filter<PointT>::PointCloud PointCloud;
	typedef typename PointCloud::Ptr PointCloudPtr;
	typedef typename PointCloud::ConstPtr PointCloudConstPtr;

public:
	/*************************************************************************/
	/**
	* \brief		ASF 构造函数 
	* \param[in]	max_window 
	* \param[in]	min_window 
	* \returns		
	* \remark		 
	*/ 
	/*************************************************************************/
	AdaptiveSurfaceFilter(double max_window = 30.0, double min_window = 1.0);
	~AdaptiveSurfaceFilter(void);

	void setParameters(const ASFParameter& param);

	/*************************************************************************/
	/**
	* \brief		获得地面点点云的序号 
	*/ 
	/*************************************************************************/
	vector<int> getGroundIndex() const { return ground_index_; }
	
	/*************************************************************************/
	/**
	* \brief		获取非地面点点云的序号 
	*/ 
	/*************************************************************************/
	vector<int> getNongroundIndex() const { return nonground_index_; }
	
	/*************************************************************************/
	/**
	* \brief	将点云分成low outlier, ground nonground	
	*/
	/*************************************************************************/
	vector<uint8_t> getClassification();

	/*************************************************************************/
	/**
	* \brief		获得滤波后的最精细一层的DEM，如果想获得质量更好的DEM，建议还是
	* 从ground point重新内插，这里获取的只是一个为了滤波临时的DEM
	* \param[out]	dem 	 
	*/ 
	/*************************************************************************/
	void getDEM(std::shared_ptr<RasterSurface>& dem) { dem = dem_; };
	/*************************************************************************/
	/**
	* \brief		如果有已知值，可以计算混淆矩阵 
	* \param[in]	ref_be reference 地面点序号
	* \param[in]	ref_obj reference 非地面点序号
	* \param[in]	flt_be 滤波的地面点序号
	* \param[in]	flt_obj 滤波的非地面点序号
	* \returns		vector<int> 混淆矩阵中的四个值，见Hu et al. 2014的表1
	* \remark		 
	*/ 
	/*************************************************************************/
	vector<int> evaluate(vector<int>& ref_be, vector<int>& ref_obj,
		vector<int>& flt_be, vector<int>& flt_obj);

private:
	double step_size_;
	double threshold_;
	double max_window_;
	double min_window_;

	double max_scale_gain_;
	double max_bending_gain_;
	double max_lambda_;

	int bending_level_;

	int m_;
	int k_;

	bool detect_low_outlier_;

	///\brief low outlier detection的参数
	///\sa LowOutlierFilter
	double low_outlier_mophological_window_;
	double low_outlier_threshold_;
	double low_outlier_mophological_size_;
	int max_low_grids_;

	vector<int> ground_index_;
	vector<int> nonground_index_;
	vector<int> low_outlier_index_;

	std::shared_ptr<RasterSurface> dem_;
protected:
	/*************************************************************************/
	/**
	* \brief		采用分段线性内插，根据bending energy计算bending gain 
	* \param[in,out]	bending 输入为bending energy的值，输出的结果是bending gain
	*/ 
	/*************************************************************************/
	void transformBendingEnergy(RasterSurface* bending);

protected:
	void applyFilter(PointCloud &output);
};
}//end oil


#ifdef PCL_NO_PRECOMPILE
#include "adaptive_surface_filter.hpp"
#endif

